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ABSTRACT 



Context. Many thermally emitting, isolated neutron stars have magnetic fields that are larger than 10'^ G. A realistic cooling model 
that includes the presence of high magnetic fields should be reconsidered. 

Aims. We investigate the effects of an anisotropic temperature distribution and Joule heating on the cooling of magnetized neutron 
stars. 

Methods. The 2D heat transfer equation with anisotropic thermal conductivity tensor and including all relevant neutrino emission 
processes is solved for realistic models of the neutron star interior and crust. 

Results. The presence of the magnetic field affects significantly the thermal surface distribution and the cooling history during both, 
the early neutrino cooling era and the late photon cooling era. 

Conclusions. There is a large effect of Joule heating on the thermal evolution of strongly magnetized neutron stars. Both magnetic 
fields and Joule heating play an important role in keeping magnetars warm for a long time. Moreover, this effect is important for 
intermediate field neutron stars and should be considered in radio-quiet isolated neutron stars or high magnetic field radio-pulsars. 

Key words. Stars: neutron - Stars: magnetic fields - Radiation mechanisms: thermal 



1. Introduction 

Observation of thermal emission from neutron stars (NSs) can 
provide not only information on the physical properties such as 
the magnetic field, temperature, and chemical composition of the 
regions where this radiation is produced but also information on 
the properties o f matter at higher densities deeper inside the star 
dYakovlev & Pe thick 2004: P age et al. 2006) . To derive this in- 
formation, we need to calculate the structure and evolution of the 
star, and compare the theoretical model with the observational 
data. Most previous studies assumed a spherically symmetric 
temperature distribution. However, there is increasing evidence 
that this is not the case for most nearby neutron stars whose ther- 
mal emiss ion is visible i n the X-ray b and of the electromagnetic 
spectrum (IZavlinll2007t lHaber]||2007l) . The anisotropic temper- 
ature distribution may be produced not only in the low density 
regions where the spectrum is formed and preliminary investiga- 
tions had focused their attention, but also in intermediate density 
regions, such as the solid crust, where a complicated magnetic 
field geometry could cause a coupled magneto-thermal evolu- 
tion. In some extreme cases, this anisotropy may even be present 
in the poorly known interior, where neutrino processes are re- 
sponsible for the energy removal. 

The observational fact that most thermally e mitting isolate d 
NSs have magnetic fields larger than 10'^ G (Haberl 2007), 
which is sometimes confirmed by spin down measurements, 
leads to the conclusion that a realistic cooling model must not 
avoid the inclusion of the effects produced by the presence of 
high magnetic fields. The transport processes that occur in the 
interior are affected by these strong magnetic fields and their 
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effects are expected to have observable consequences, in par- 
ticular for highly magnetized NSs or magnetars. Moreover, the 
large surface magnetic field strengths inferred from the ob- 
servations probably indicate that the interior field could reach 
even larger values, as theor etically predicted by some models 
([Thompson & Duncanlll993h . 

The presence of a magnetic field affects the transport proper- 
ties of all plasma components, especially the electrons. In gen- 
eral, the motion of free electrons perpendicular to the magnetic 
field is quantized in Landau levels, and the thermal and electri- 
cal conductivities exhibit quantum oscillations. In the limit of 
a strongly quantizing field, in which almost all electrons popu- 
late the lowest level, such as in the envelope of a NS, a quan- 
tum description is necessary to calculate th e thermal and elec- 
trical c ond uctiv i ties. E arlier calculations by ICanuto & Chiuderil 
('1970') and lltohl (Il975h concluded that the electron thermal con- 
ductivity is strongly suppressed in the direction perpendicular to 
the magnetic field and increased along the magnetic field lines, 
which reduces the thermal insulation of the envelope [heat blan- 
keting). Thus, there is an anisotropic heat transport in the NS's 
envelope governed by the magnetic field geometry, that produces 
a non-uniform surface temperature. 

The anisotropy in the surface temperature of a NS ap- 
pears to be confirmed by the analysis of observational data 
from isolated NSs (see IZavlinl ("2007) and Haberl (2007) for 
reviews on the current status of theory and observations). 
The mismatch between the extrapolation to low energy of 
fits to the X-ray spectra, and the observed Rayleigh-Jeans 
tail in the optical band (optical excess flux), cannot be ad- 
dressed using a unique temperature. Several simult aneous fits 
to rn u ltiwavelength spectra of RX J 1 85 6.5-3754 (Poiis et all 
I2OO2I: iTriimper et al] l2004l) . RBS 1223 dSchwope et alJ I2005L 
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120071) . and RX J0720.4-3125 dPerez-Azorin et al.1 l2006bl) are 
explained by a small hot emitting area ^ 10-20 km^, and an 
extended cooler component. Another piece of evidence that 
strongly supports the nonuniform temperature distribution are 
pulsations in the X-ray signal of some objects of amplitudes ^ 
5-30 %, some of which have irregular light curves that point to- 
wards a non-dipolar temperature distribution. All of these facts 
reveal that the idealized picture of a NS with a dipolar magnetic 
field and uniform surface temperature is oversimplified . 

In a pioneering work, ' Greenstein & Hartk^ ( 1198 3^ obtained 
the temperature at the surface of a NS as a function of the mag- 
netic field inclination angle in a simplified plane-parallel approx- 
imation. This model was applied to diff'erent magnetic field con- 
figurations and the observational consequences of a non-uniform 
temperature distribution were analyzed in the pulsars Vela 
and G eminga among others (Page 19951). [Potekhin & YakovlevI 
( 1200 lb improved the former calculations including realistic ther- 
mal conductivities. Nevertheless, the temperature anisotropy as 
generated in the envelope may be insufficiently to be consis- 
tent with the observed thermal distribution and, in this case, 
should originate deeper inside the NS (Geppert et al. 200^ 
iPerez-Azorin et al.ll2006al) . 

Crustal confined magnetic fields could be responsible for 
the surface thermal anisotropy. In the crust, even if a strong 
magnetic field is present, the electrons occupy a large num- 
ber of Landau levels and the classical approximation remains 
valid during a long time in the thermal evolution. The mag- 
netic field limits the movement of electrons in the direction 
perpendicular to the field and, since they are the main car- 
riers of the heat transport, the thermal conductivity in this 
direction is highly suppressed, while remaining almost unaf- 
fected along the field lines. Temperature distributions in the 
crust were obtained as stationary solut ions of the diffusion 
equation with axial symmetry (Geppert et al.l l2004l) . The ap- 
proach assumes an isothermal core and a magnetized envelope 
as an inner and outer boundary condition, respectively. The re- 
sults show important deviations from the crust isothermal case 
for crustal confined magnetic fields with strengths larger than 
10'^ G and temperatures below 10** K. Similar conclusions 
were obtained considering not only poloidal but also toroidal 
components for the magnetic field dPerez-Azorm et akl l2006at 
iGeppert et alj|2006l) . This models succeeded in explaining si- 
multaneously the observed X-ray spectrum, the optical excess, 
the pulsed fraction, and other sp ectral features for some iso lated 
NS such as RX J0 720.4-3125 dPerez-Azorm et alj |2006bh and 
RX J1856.5-3754 (iGeppert et alJl2006l) . 

Non-uniform surface temperature in NSs was studied by dif- 
fere nt authors using simplified models (Shib anov & YakovlevI 
11996; Po tekhin & Yakovlev 2001). Although these models can 
provide useful information, a detailed investigation of heat trans- 
port in 2D must be completed to obtain more solid conclu- 
sions. However, this is not the only effect that must be re- 
visited to study the cooling of NSs. For isolated NSs, differ- 
ent r elevant magnetic field dissipatio n processes were identi- 
fied dGoldreich & Reiseneggeii Il992h . The Ohmic dissipation 
rate is determined by the finite conductivity of the constituent 
matter In the crust, the electrical resistivity is mainly due 
to electron-phonon an d electron-impurity scattering processes 
dFlowers & Itoh|[r976l) . resulting in more efficient Ohmic dissi- 
pation than in the fluid interior. The strong temperature depen- 
dence of the resistivity leads to rapid dissipation of the magnetic 
energy in the outermost low-density regions during the early 
evolution of a hot NS, which becomes less relevant as the star 
cools down. Joule heating in the crustal layers due to Ohmic de- 



cay was thought to affect only the late photon cooUng era in old 
NS (> 10^ yr), and to be an efficient mechanism to maintain 
the surface temperature as high as ^ lO'*"^ K for a long time 
(iMiralles et al.lll998l) . [Pa"ge et alldlOOO ) studied the 1-D thermal 
evolution of NSs combined with an evolving Stokes function that 
defines a purely poloidal, dipolar magnetic field. The Joule heat- 
ing rate was evaluated averaging the currents over the azimuthal 
angle. However, for strongly magnetized NS, Joule heating can 
be important much ea rlier in the evolution. In a recent work, 
iKaminker et al.l (12006^ placed a heat source inside the outer crust 
of a young, warm, magnetar of field strength 5 x 10'^ G. To ex- 
plain observations, they concluded that the heat source should 
be located at a density < 5 x 10" g cm"-', and the heating rate 
should be > 10^** erg cm"^ s"' for at least 5 X 10"* yr. Anisotropic 
heat transport is neglected in these simulations, which were per- 
formed in spherical symmetry, assuming that it will not affect the 
results in the early evolution. Nevertheless we will show that, in 
2D simulations, the effect of anisotropic heat transport is impor- 
tant. 

In addition to purely Ohmic dissipation, strongly magnetized 
NSs can also experience a Hall drift with a drift velocity pro- 
portional to the magnetic field strength. Although the Hall drift 
conserves the magnetic energy and it is not a dissipative mecha- 
nism by itself, it can enhance the Ohmic decay by compressing 
the field into small scales, or by displacing currents to regions 
with higher resistivity, where Ohmic dissipation is more effi- 
cient. Recently, the first 2D-long term simulations of the mag- 
netic field evolution in the crust stud ied the interplay of Oh mic 
dissipation and the Hall drift effect (iPons & Geppert! 120071) . It 
was shown that, for magnetar field strength, the characteristic 
timescale during which Hall drift influences Ohmic dissipation 
is of about lO'* yr. All of these studies imply that both field decay 
and Joule heating play a role in the cooling of neutron stars born 
with field sti-engths > lO'^ G. 

We will show that, during the neutrino cooling era and the 
early stages of the photon cooling era, the thermal evolution is 
coupled to the magnetic field decay, since both cooling and mag- 
netic field diffusion proceed on a similar timescale (» 10^ yr). 
The energy released by magnetic field decay in the crust could 
be an important heat source that modifies or even controls the 
thermal evo lution of a NS. O bservational evidence of this fact 
is shown in iPons et al.l d2007l) . They found a strong correlation 
between the inferred magnetic field and the surface temperature 
for a wide range of magnetic fields: from magnetars (> 10'^ 
G), through radio-quiet isolated neutron stars (^ 10'^ G) down 
to some ordinary pulsars (< 10'^ G). The main conclusion is 
that, rather independently from the stellar structure and the mat- 
ter composition, the correlation can be explained by the decay of 
currents on a timescale of ^ 10^ yr 

The aim of the present work is to study in a more consis- 
tent way the cooling of a realistic NS under the effects of large 
magnetic fields, including the effects of an anisotropic tempera- 
ture distribution and Joule heating in 2D simulations. As a first 
step towards a fully coupled magneto-thermal evolution, a phe- 
nomenological law for the magnetic field decay is considered. 

This article is structured as follows. In Sect. 2 we discuss 
the equations governing the magnetic field structure and evo- 
lution, while Sect. 3 is devoted to the thermal evolution equa- 
tions. Sect. 4 presents the microphysics inputs. Sect. 5 and 6 
contain our results for weakly and strongly magnetized NSs, re- 
spectively. In Sect. 7, we focus on the effects of field decay and 
Joule heating on the cooling history of a NS. Finally, in Sect. 8 
we present the main conclusions and perspectives of the present 
work. 
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2. Magnetic field structure and evolution. 



2. 1 . Magnetic field geometry 



While the large-scale external structure of the magnetic field of 
NSs is usually represented by the vacuum solution of an external 
dipole, or sometimes a more complex magnetosphere, the struc- 
ture of the magnetic field in the interior of NSs is poorly known. 
Results from MHD simulations show that stable configurations 
require the coexistence of both poloidal and toroidal compo- 
nents, approximately of the same strength (Braithw aite & Spruill 
[2004), although predominantl y poloidal configurations ma y be 
stabilized by rapid rotation jGeppert & Rheinhardtl l2006h . In 
general, a realistic NS magnetic field model should contain both 
components, and their location and relative strength should vary . 
Moreover, two-dimensional simulations (Pons & Geppert 2007") 
showed that, while the initial magnetic field configuration deter- 
mines the early evolution of the field {t < 10^ yr), at later stages a 
more stable configuration, consisting of a dipolar poloidal com- 
ponent and a higher order toroidal component, is preferred. 

We consider the Newtonian approximation of a NS magnetic 
field because general relativistic corrections are not important in 
our study. In axial symmetry, the magnetic fi eld can be dec om- 
posed into poloidal and toroidal components (lRaedleij|2000h 



B = B 



pol 



(1) 



which are represented, respectively, by two scalar functions S, 
t: 



Bpoi = Vx(rxV^) 



(2) 
(3) 



Here, S and T" depend on the spherical coordinates r, 0, and r is 
a radial vector 

Expanding the angular part of the scalar functions in 
Legendre polynomials, we can write 

P,(cos 6) 



S(r, e)^C 2^ S,{r, t) , 



- P/(cos0) 



where P;(cos 6) is the Legendre polynomial of order I and C is a 
normalization constant. For / = 1, wich represents dipolar fields, 
after normalizing the field to its surface value at the magnetic 
pole, B, (C - R^^ B/2) and the radial coordinate to the NS radius 
(x = r/R^s), the components of the magnetic field can be written 
in terms of the two functions, Si{x, t) and T\{x, t), as follows 



Br = B- 



COS0 



Be 

Bd, 



-B 



-Si{x,t) 

x~ 

sinO dS](x, f) 



2x 



dx 

sin0„ 
B—Ti(xj), 
Lx 



(5) 



where in the following we omit the subindex (I - V) for clar- 
ity. We note that S(x, t) is normalized such that it reaches the 
value of 1 at the surface. These two arbitrary functions are sub- 
ject to suitable boundary conditions. To match continuously the 
external vacuum dipole solution, for example, S{x, t) must sat- 
isfy dS{x, t)ldx - -S{x, t), at X = 1. Other boundary conditions 
are discussed below. 



From the above general form of the magnetic field components, 
different interesting cases can easily be recovered. We describe 
three possible configurations that we explored in this work. 

2.1 .1 . Force-free fields (FF model) 

One of the particular models we consider here are the force-free 
fields. They satisfy: 



VxB ^ fiB, B -V/j^O 



(6) 



where yu is a parameter related to the magnetic field curva- 
ture, which naively can be interpreted as a wavenumber of the 
Stokes function S. For simplicity, we consider solutions with 
/i =constant such that the second equation is satisfied automati- 
cally. A general interior solution that fulfils the equality between 
the two components (r, 6) in the first equation can be obtained 
by choosing 



T(x, f) = ^S{x, t). 



(7) 



Factoring the time dependence in an arbitrary function, S{x, t) - 
f(t)A{x), the 0-component of the first equality in Eq.|6]produces 
a form of the Riccati-Bessel equation for A(x) whose solution 
can be written analytically in terms of the spherical Bessel func- 
tions of the first and second kind. For I - 1, we have 

A{x) — axji{x) + bxni{x) , 
sin X cos X 



ji(x) = ,2 
ni(x) = - 



x^ 

cos X 



X 

sinx 



X^ X 

where x - fiRf^^x. From this, the magnetic field is given by 

cosO 
B, = B^-A{x), 



(8) 



B„ = -B 



smd dA{x) 



2x dx 
sinO 

(4) = BfiRj^s-^A(x) . 



(9) 



This family of solutions is parameterized by B and the value of 
the dimensionless quantity /i7?NS- To match continuously the ex- 
ternal vacuum dipole solution, one must choose a - cos(yL(/?Ns), 
b = sin(;u/?Ns)- 

If the magnetic field extends to the center of the NS, only 
the regular solutions at x = (ji) must be considered, i.e., we 
must set b = 0, which directly determines fi. Due to the su- 
perconducting nature of the fluid core, the magne tic field may 
be expelled and confined to the crustal region (iJonesll 19871: 
iKonenkov & Gepperlll2001h . 

This is of course a simplification, since in a type II supercon- 
ductor the magnetic field would be organized in flux tubes with 
complex geometries, but it suffices for our purposes to establish 
qualitative differences between core- and crustal-fields. In the 
case of magnetic fields confined to the crustal region, from the 
core radius (Rcok) to 7?ns, one must adjust jj to have a vanishing 
radial component in the crust-core interface. This can be done 
by solving 



tan [yU (/^core - ^NS)] = /"^c 



(10) 



The values of the parameter fi obtained for the NS models used 
in this paper are listed in Table 14.1] In Fig.[T]we show the three 
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FF 

-- TCI \ 

TC2 \ 




Fig. 1. Normalized magnetic field components for the force free 
case in the crust: Br/B/ cos9 (dashed line), Bg/B/ sinO (dashed 
dotted), and B^/B/ sin (solid line) vs normalized radial coordi- 
nate X. 



Fig. 2. Normalized toroidal field components B^/B/ sinO vs 
normalized radial coordinate x in the crust. Three different mod- 
els are shown: FF (solid lines), TCI (dashed lines) and TC2 
(dashed dotted lines). 



normalized components of the crustal confined force free field 
for the LM model. 

This force-free solution can easily be extended to higher or- 
der multipoles, e.g. quadrupole, by replacing the angular depen- 
dence by the corresponding Legendre polynomial and using the 
corresponding spherical Bessel functions of the same index I. 
From the above general expression of force-free fields, some of 
the cases usually considered in the literature can be recovered. 

2.1 .2. Configurations witli oilier toroidal fields (TC1 and TC2) 

We consider another two models with crustal-confined toroidal 
fields that obey 

Tixj) = TQx(l-xf(x-R,,,JRt,s) (Model TCI) (11) 
T(x, t) = Tqx{1 - x){x - R,orJRNsf' (Model TC2) , (12) 

with the same poloidal component as in the FF case. The con- 
stant T() is fixed such that B^ is one order of magnitude larger 
than B,- at the NS surface. In the latter two configurations, the 
maximum of B^ is close to the crust-core boundary or close to 
the crust-envelope boundary, respectively. The resulting radial 
profiles of B^ are shown in Fig. |2] We note that the toroidal 
component of the FF configuration penetrates into the envelope, 
while the remaining two (TCI and TC2) are confined to the 
crust. 

2.1 .3. Crustal poloidal fields (PC model) 

If we assume that the magnetic field is confined to the crust, 
and maintained by purely toroidal currents, we can simply set 
T{x, f) = and, from Eq. |5] we have 



costy 
Br = B—^S(x,t), 



Bf, 



-B 



sine dS(x,t) 
2x dx 



(13) 



where again the boundary conditions dS{x,t)/dx - -S(x,t) 
at X = 1, and S{x,t) = at x - R^ore must be fulfilled. 
In general, S(x, t) does not need to coincide with the function 
A{x) expressed above in terms of the spherical Bessel functions. 
However, given the freedom in the choice of S{x, t), we prefer 
to use the analytical form of A(x) to determine the poloidal field. 



rather than specifying a similar solution that would be equally 
arbitrary. 

2.1.4. Core dipolar solutions (CD model) 

The extension of the vacuum solution towards the interior can 
be shown to correspond to the limit yU — > of the non-regular 
function n i , explicitly. 



Br = B- 



cos 6 



Bo = -B 



, sin0 
2x^' 



Bs^O. 



(14) 



Although this solution diverges at jc = 0, it has been used in the 
literature to represent the magnetic field structure in the crust, 
assuming that the field is reorganized in an unknown form in the 
core. Alternatively, one can also take the limit yU — > of the reg- 
ular spherical Bessel function ji. This leads to a homogeneous 
field aligned with the magnetic axis. 

2.2. Field decay and Joule heating 

The induction equation that describes the evolution of the mag- 
netic field in the crust is 



dB 

5f 



— = -Vx 



r]V xB + 



Anen, 



(VxB)xB 



(15) 



where 77 = ^ is the electrical resistivity, cr is the electrical 
conductivity parallel to the field lines, is the electron density, 
and e the electron charge. 

The first term in the bracket is purely diffusive (Ohmic) and 
the second corresponds to the Hall term. Taking the scalar prod- 
uct of B by Eq. ( fTSl l, and integrating over the volume, it can be 
seen that the Hall term does not contribute to the dissipation of 
energy, but it redistributes the magnetic energy from one place to 
another A force-free field satisfying VxB - fiB is not subject to 
the Hall term and, if 77 is constant throughout the crust volume, 
the induction equation is reduced to 



(16) 



which shows that purely Ohmic dissipation is exponential and 
proceeds on a typical timescale rohm = ^/w^- 
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In a realistic case the situation is more complex, since the 
non-linear evolution of the Hall term must be taken into account 
and the conductivity and electron density profiles are not con- 
stant. Even if we start from a force-free magnetic field, the effect 
of a resistivity gradient leads to a fast modification of its geom- 
etry, and the Hall term becomes immediately important. 

For simplicity, and with the purpose of investigating qual- 
itatively the effects of magnetic field decay, we include phe- 
nomenologically a first stage with rapid (non-exponential) de- 
cay, and a late stage with purely Ohmic dissipation (exponen- 
tial). We assume that the geometry of the field is fixed and the 
temporal dependence is included only in the normaUzation value 
B according to 



B^Bo 



exp(-f/Tohm) 



1 + ^(1 



•exp(-f/Tohm)) 



(17) 



where rohm is the Ohmic characteristic time, and the typical 
timescale of the fast, initial stage is defined by THaii- This is the 
analytical solution of the differential equation 



dB 

dt 



B 



1 B^ 



TOhm BqT Hall 



(18) 



that takes into account the approximate dependence of the Hall 
timescale on the magnetic field ^/B^). We note that THaii 
should be interpreted as the Hall timescale corresponding to the 
initial magnetic field strength Bo- In the early evolution, when 

t Tohm, 



B-Bo(l+f/THall)"' 

while for late stages, when / > rohn 
B ^ Bo exp(-f /rohm) 



(19) 



(20) 



This simple law reprod uces qualitatively the r esults from 
more complex simulations jPons & GeppertI [20071) and facili- 
tates the implementation of field decay in the cooling process 
of NSs for different Ohmic and Hall timescales, treated as sim- 
ple constant parameters. The initial Hall stage, in which the Hall 
drift qualitatively affects the thermal evolution, is of particular 
importance for models of highly magnetized NS, e.g, for mag- 
netars the field can dissipate 75% of the energy in x tum- In 
contrast, the late Ohmic stage lasts for about Tohm - 10^ yr. 

If the field is anchored into the superconducting core, the 
results will be different. It is not the purpose of this paper to dis- 
cuss such a possibility, which deserves a separate analysis, but 
to investigate the possible effects of crustal fields that enhance 
the surface temperature anisotropy and are subject to Ohmic dis- 
sipation and, consequently. Joule heating. 



3. Thermal evolution 

3.1. The diffusion equation in axial symmetry 

Assuming that deformations with respect to the spherically- 
symmetric case due to rotation, magnetic field, and temperature 
distribution do not aff ect the metric in th e interior of a NS, we 
use the standard form dMisner et al 

ds^ = -e^'^dt^ + e^^dr^ + p-dQ} . (21) 

Using this background metric but considering an axially- 
symmetric temperature distribution, the thermal evolution of a 
NS can be described by the energy balance equation 



Cve"^+V-(e2'i'F) = e--2 
ot 



„2<t. 



(22) 



where c,. is the specific heat per unit volume and Q is the energy 
loss/gain by v-emission. Joule heating, accretion heating, etc.. In 
the diffusion limit, the heat flux is simply 

F = -e-'^k ■ V(e*r) (23) 

where k is the thermal conductivity tensor Defining the red- 
shifted temperature to be T = e^T, the components of the red- 
shifted flux F = e^'^F can be written explicitly as follows 



F,^-e''[K,,e-^d,-T+'-fdeT] 



(24) 



where the 0-component is not relevant because of the axial sym- 
metry. 

The total conductivity tensor, k, must include the contribu- 
tions of all relevant carriers, which are of interest in the solid 
crust: electrons, neutrons, protons and phonons 

k^ke + k„ + kp + kph . (25) 

The heat is transported primarily by electrons, which provide 
the dominant contribution. Radiative transport is important close 
to the surface, but the outer region is considered by means of 
boundary conditions (Sect. [372] l in place of direct calculation. 

For magnetized NS, the electron thermal conductivity tensor 
becomes anisotropic: in the direction perpendicular to the mag- 
netic field, its strength is strongly diminished, which causes a 
suppression of the heat flow orthogonal to the magnetic field 
lines. The ratio of conductivities along and orthogonal to the 
magnetic field can be defined in terms of the magnetization pa- 
rameter, (Obt, as 



(26) 



where r is the electron relaxation time ('Urpin & YakovlevI 
1980), and a>B is the classical electron gyrofrequency corre- 
sponding to a magnetic field strength B 

eB 

(^B^ — , (27) 

where m* is the electron effective mass. The dimensionless quan- 
tity (ji)bt is an indicator of the suppression of the thermal conduc- 
tivity in the transverse direction. When cjbt » 1, the effects of 
the magnetic field on the transport properties are crucial. 

In spherical coordinates, and choosing the polar axis to coin- 
cide with the axis of symmetry of the magnetic field, the electron 
contribution can be written as follows 



I + {ojbtY 



( h 

bre 

brd> 



bre 
bee 
bffd, 



br^ 

be^, 

bdidi 



■ (^bt 





-b^ 

bf, 





-br 



-be 
br 




(28) 

where / is the identity matrix, and b,-, bg, b^ are the components 
of the unit vector in the direction of the magnetic field, and 
bij - bjbj for i, j = r, 6, <p. Using the above expression for kg, 
the electron part of the flux reads, in closed form: 

F, = -e'^K^ [Vf -I- (wbt)2 {b-VT)-b+ cobt [b x Vf )] . (29) 

The thermal evolution Eq. (l22T i. with the above expression 
for the fluxes, is solved numerically for a given background mag- 
netic field with fixed geometry and strength that varies with time 
according to Eq. ( [TtI i. The emissivity terms on the right-hand 
side of Eq. (l22T i and the specific heat of the first term of the same 
equation are considered in the next section. 
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3.2. Boundary conditions 

For numerical reasons, the thermal evolution equation is diffi- 
cult to solve in the thin layer which consist of the envelope, of 
a few meters depth, and the atmosphere, of a few centimeters, 
in which radiative equilibrium is established and the observed 
spectrum is generated. Since this outer layer has a small scale 
and its thermal relaxation time is much shorter than the overall 
evolutionary time, the usual approach is to use results of station- 
ary, plane-parallel, envelope models to obtain a phenomenolog- 
ical fit that relates the temperature at the bottom of the envelope 
Tb, with the surface temperature Tj. This ''Th-Ts relationship" 
can be used to implement boundary conditions, because the sur- 
face flux can then be calculated for a given temperature at the 
base of the envelope, which corresponds to the outer point of the 
numerical grid in our cooling simulations. 

Models assuming a non-magnetized envelope made of iron 
and iron-like nu clei show that the surface te mperature is related 
to Th as follows jGudmundsson et aTll 19831) 



n,8 = 1.288 



7^4 

s,6 



;i4 



0.455 



(30) 



where gi4 is the surface gravity in units of 10''* cm s ^, 7/,^^, is 
Th in 10** K, and 7^,6 is in 10^ K. 

Since the magnetic field increases the heat permeability of 
the envelope in regions where the magnetic field lines are radial 
but strongly suppresse s it where the rnagne tic field lines are al- 
most tangential (Potek hin & Yakovlevl200lh . this implies a large 
anisotropic distribution of T^, which depends on the magnetic 
field geometry. In iron magnetized envelopes, the surface tem- 
perature depends on the angle (f that the magnetic field forms 
with respect to the normal to the NS surface by means of a func- 
tion X: 



T,(B, if, g, Th) * rf (g, Th) X{B, if, Th) 



where 



(31) 



(32) 



and C = 0.iri,8 -0.001 g^^^ y/oTT^. The function was fitted 
by decomposing into transversal and longitudinal parts as 



X(B, (fi, Th) = [ Xl/-(B, Th) cos^ ip + X^^^iB, Th) sin^ (pf^ , 

which is vahd for B < 10'^ G and 10^ K < T,, < lO**^ K, with 
the additional constraint that T^ > 2 x 10^ K. 

Strongly magnet ized envelopes were revisited by 
iPotekhin et al.l (l2007h . who reconsidered neutrino emission 
processes that are activated by strong fields, that is neutrino 
synchrotron. These processes were found to lower the surface 
temperature at a fixed Th. To take this effect into account, we 
introduce a maximum surface temperature Tf^'^((p) that can be 
reached for a given Th, which we parameterize as a function of 
B to reproduce their results. 

In this work, we assume an iron composition for the envelope 
and focus on the magnetic corrections to the transport due to the 
presence of a large field. Nevertheless, if light elements were 
present, which is very unlikely because the large magnetic field 
suppresses accretion, they strongly reduce the blanketing effect 
and the relations used here should be revised. Another possibility 
is that the gaseous atmosphere and the outer envelope conden- 
sates to a solid state due to the cohesive interaction between ions 



Table 1. Central density pc, mass M, radius Rf^s, crust thickness 
A/?crust> and jj parameter for the two cases used in this work: low 
mass (LM) and high mass (HM). 



Model 


Pc 


M 










(g cm~') 




(km) 


(km) 


(km-') 


LM 


8.1 10" 


1.35 


12.83 


1.24 


1.32 


HM 


1.1 10'^ 


1.63 


12.36 


0.86 


1.87 



caused by the magnetic field. This condensed surface has differ- 
ent emission properties and, consequently, t he boundary condi- 
tion mu st be recalculated, as for example in iPerez-Azorm et all 
(l2006ah . This scenario will be studied in future work. 



4. Microphysics inputs 

4.1. EoS and superfluidity 

To build the background NS model, we used a Skyrme-type 
equation of state (EoS), at zero temperature to describe both the 
NS crust and the liquid core, based on t he effective nuclear in- 
teraction SLy (Douchin & Haens"ell2001h . The low density EoS, 
below the neutron drip point, employed is that of Baym et al. 
(1971). Throughout this work we use two models: a low mass 
neutron star (LM) and a high mass (HM) NS, which have the 
properties listed in Table |4.1] For the chosen EoS, the crust-core 
interface is at 0.46 po, where po - 2.8x 10'^ g cm""* is the nuclear 
saturation density, and for both models the crust thickness is ap- 
proximately 1 km, which defines a characteristic length scale for 
the confinement of the crustal magnetic field. 

In Fig. |3] the number of particles per baryon {Y(„^p^i,)) and 
the fraction of nucleons inside heavy nuclei (X/,) are shown as 
a function of the density. In the upper horizontal axis, the scale 
shows the value of the radial coordinate that limits each region: 
the outer and inner crust, and the outer and inner core, for the LM 
and HM NS models. We do not include muons in our equation 
of state. 

Pairing in nuclear matter can play an important role in NS 
cooling, without affecting significantly the EoS, but strongly 
modifying neutrino emissivities and specific heat. In fact, for 
the paired component, these are suppressed by exponential 
Boltzmann factors. If superfluidity (SF) occurs inside NSs, i.e. 
when T is below a critical temperature (r^), the inclusion of 
these suppression factors will have important consequences on 
the thermal evolution as we see in the following subsections. 
We consider the pairing of neutrons in the crust in the n ^So 
state, protons in the core in the p ^ Sp state, and the « s tate, 
for neutrons in the cor e. Following iKaminker et al.l (1200 lb and 
lAndersson et al.l (l2005h . we use a phenomenological formula for 
the momentum dependence of the energy gap at zero tempera- 
ture 



A(kF,N) = Ao 



(kF,N - ^2) 



(kf^N - ko)^ + k] (kF,N - k2)^ + H 



(33) 



where kF,N = (37r rti^) ' is the Fermi momentum and is the 
particle density of each type of nucleons (A^ = n,p) involved. 
The parameters Ay and ki, i = 1..4 are values fitted to recent 
model calculations listed in Table |2] This expression is valid for 
^0 < kp^N < k2, with vanishing A outside this range. The density 
dependence of the gaps is plotted in Fig.|4] 

For the n 'S'o pairing, t he bare interaction predicts a max- 
imum gap A™"'" ^ 3 MeV dSchulze et al.lll998b . but the polar- 



ization effects reduce it by a factor 2-3, giving A" 



1 MeV 
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Fig. 3. NS composition for the EoS employed. Number of parti- 
cles per baryon as a function of the density: (solid line), Yp 
(solid line with plus symbols), and F,, (dashed line). Xj, is indi- 
cated by dashed-dotted lines. The scale in the upper horizontal 
axis indicates the corresponding radial coordinate at each den- 
sity for both LM and HM models. 



Table 2. Parameterization and references of the energy gaps for 
superfluid states 



Label 


Ao 


^0 


h 




h 


Ref. 




(MeV) 


(fm-') 


(fm-') 


(fm-') 


(fm-') 




n 'So 


a 


68 


0.1 


4 


1.7 


4 


1 


b 


4 


0.4 


1.5 


1.65 


0.05 


2 


c 


22 


0.3 


0.09 


1.05 


4 


3 


p'So 


e 


61 





6 


1.1 


0.6 


4 


f 


55 


0.15 


4 


1.27 


4 


5 


n^P2 


h 


4.8 


1.07 


1.8 


3.2 


2 


6 


k 


0.42 


1.1 


0.5 


2.7 


0.5 


7 


m 


2.9 


1.21 


0.5 


1.62 


0.5 


4 



References. (1) IWambach et alj jl993h: (2) 
([1998); (3) IChenetal.1 ( Il986h- 
(5) Amundsen & Ost gaardI d 19851) 
(7) lElgar0vet al.n l 996bl) 



(4) lElgar0v et ^ 
(6) llSdoetS 



ISchulze et alj 
' ( Il996j) : 
(Hsi); 



at kf n ^ 0.7 - 0.8 fm-'. This is the case for the approximation 
a, which in our NS models peaks at p = 4 x lO'"* g cm--' in 
the inner crust (Fig.|4li. Although model calculations have found 
some agreement about the value of the maximum energy gaps, 
its precise location is uncertain and may vary in the different ap- 
proaches, as shown in cases b and c. The corresponding critical 
temperatures for the .s-wave can be calculated approximately to 
be Tc — 0.56 A(r = 0), which implies a maximum for neutrons 
of T^^ = 9 X 10'' K, for the case a. As shown later, this high 
temperature implies that neutrons become superfluid in the crust 
during the early cooling of a NS and the most important conse- 
quence is that the crustal specific heat is suppressed. 

The calculations for the p '^o pairing take into account 
the presence of the neutron gas and depend also on the proton 
fraction through the symmetry energy of the EoS. For differ- 
ent approaches, such as case e and /, A™^" is located at about 




le+15 



p (g cm ^ 



Fig. 4. Energy gaps for superfluidity as a function of the density. 
Lines denote: case a (thick solid), case b (short dashed), and 
case c (long dashed) for n'5o; case e (dashed dotted) and case / 
(double dashed dotted) for p'^So; case h (thin solid) and case k 
(dotted) for «^P2- The right axes show Tc for 'S'o and pairing 
states. 



kf p ^ 0.4 - 0.5 fm-', which is much smaller than for neutrons, 
due to the smaller proton effective mass. Nevertheless, due to the 
proton number density, the peak is shifted top-2xl0''*g cm-^ 
in the outer core of our NSs, as shown in Fig. H) Most models 
agree that the proton energy gap should vanish at kf p > 1.5 
fm-', i.e. at high densities inside the star p > 10'^ g cm--'. For 
the cases considered here, T™^" ^ 2-6 x 10^ K, indicating that 
also proton superfluidity will be present from the very beginning 
of the NS thermal evolution. Due to the charge of the protons, 
the superfluid is also in a superconducting (SC) state. 

The situation for the pairing of neutrons in the core is more 
complicated, because the state has coupled anisotropic gap 
equations. Some calculations show that the energy gap should 
be reduced by a factor of 2-3, as in the proton case, due to 
the lower neutron effective mass in very dense matter. But rel- 
ativistic effects become important and there is no conclusive ap- 
proach to the problem. Thus, in our calculations we considered 
three different cases that reflect this uncertainty: h, k and m with 
A™" ^ 0.6,0.1,0.02 MeV, respectively (Table |2l). The location 
of the maximum varies as well, at fe/7„ ^ 1.4-2fm-', i.e. atp ^ 2- 
6 X 10'^ g cm--', as plotted in Fig.|4l in which case m is omitted 
because it is not visible in this scal e. We note that for the p-wave 
we take that = 0.82 A(r = 0) (B ailin & Lov3l 19841) . which 
corresponds to a a wide range of T™" ^ 2 x 10^-6 x 10^ K for 
the chosen models. 

The temperature dependence of the e nergy gap that we use is 
the approximate functional form given bv lLevenfish & YakovlevI 
(119941) : 



MT) 
kT 



Tr 



P 



J 



^ITJTc TIT, 



(34) 



where a = 1.456,/? = 0.157, and y = 1.764 for '^o, and a = 
0.189, /3 = 0, and y = 1.188 for states. 

From these considerations, it is clear that a NS at the be- 
ginning of its cooling history should contain superfluid neutrons 
in the crust and superconducting protons in the core, while the 
occurrence of neutron pairing in the core is rather model depen- 
dent. 
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4.2. Thermal conductivity 

In NS cooling simulations, the thermal conductivity should be 
calculated over a region covering a large range of densities, from 
the core 10'^ g cm"-') to the outer crust 10^ g cm"-'). 

Schematically, the thermal conductivity tensor can be written 
for e ach carrier in terms o f the effective relaxation time tensor, 
fpff (iFlowers & Itohlll97d; lUrpin & Yakovlevlll980l: lltoh et alj 
liaSJ), as follows. 



3m* 



-Teff, 



(35) 



where n is the carrier number density, m* is its effective mass, 
and feff is a tensor whose components are interpreted as effec- 
tive relaxation times. In the non-quantizing case, these relax- 
ation times can be written in terms of the non-magnetic relax- 
ation time, which is calculated to be the inverse of the sum of all 
collision frequencies of the processes involved. 

In the inner liquid core, we include contributions from 
electrons, neutro ns and protons (Gnedin & Yakovlev 119951; 
iBaiko et aLll200ll) . without taking account of the effects of the 
magnetic field because of proton superconductivity: the field is 
either expelled from the core or confined into flux tubes that oc- 
cupy a small fraction of its volume. We note that, if the magnetic 
field does not affect transport properties, a large thermal conduc- 
tivity of matter is produced soon after birth in an isothermal core 
(Fig- 13, which implies that the precise value of the thermal con- 
ductivity is not important. 

In the solid crust, only electron and phonon transport are 
considered. While phonon conductivity is negligible in non- 
magnetic neutron stars, this situation changes when the mag- 
netization parameter becomes large. Since electron transport is 
drastically suppressed in the direction transverse to the magnetic 
field, the phonon contribution may become dominant at low den- 
sity as shown in Fig.|5j 




le+lO le+11 le+12 le+13 le+U 

P (g cm" ) 



.1 



T'ZJi I I 1 1 

le+10 le+11 le+12 le+13 le+14 

P (g cm" ) 



Fig. 5. Thermal conductivity contributions as a function of the 
density for different fixed temperatures, /cf is shown with solid 
lines, with dashed lines (thick for B - 10'^ G and thin for 
B - lO'"^ G), and Kph with dotted dashed lines. 



In our calculations, we use the non-quantizing electron con- 
ductivities from the public code of A. Potekhin (1999fl The 

' www. iof f e . rssi . ru/astro/conduct/condmag . html 



three electron scattering processes that play a role in our scenario 
are scattering off ions, electron-phonon scattering, and scatter- 
ing off impurities. Semi-analytic expressions and fitting formu- 
lae for the relaxation time and thermal conductivity along the 
magnetic field for all three processes, were derived by Potekhin 
& Yakovlev (1996). 

At high temperatures, the phonon conductivity of the lattice 
is determined mainly by Umklapp processes, and can be approx- 
imated by the expression 



1 



(36) 



where is the sound speed, c,, the specific heat, and Apf, the 
phonon mean free path in the lattice. In Fig. |5] it can be seen 
that the phonon contribution becomes more important at lower 
densities as the temperature decreases and the liquid solidifies 
into a lattice. 

Chugunov & Haensel IChugunov & Haensell (l2007h revised 
the ion thermal conductivity in neutron star envelopes. They in- 
cluded the contribution of electron-phonon scattering and im- 
proved the calculations of phonon-phonon scattering. Our esti- 
mates for Aph are larger than their results by a factor of a few, 
depending on the density, which results in a smaller temperature 
anisotropy. However their results are more applicable to the neu- 
tron star envelope, than for the crust. The main reason is that, at 
temperatures smaller than the Debye temperature, the inclusion 
of the effect of impurities and defects in the crystal becomes nec- 
essary. Given our limited knowledge of the impurity content of 
the inner crust, which may affect the results, we do not include 
phonon-impurity interactions in our simulations. In principle, its 
effect would be to reduce the phonon mean free path, but it is 
unclear how to calculate accurately this contribution at low tem- 
perature. 

4.3. Specific iieat 

In normal non-superfluid neutron star matter, most of the to- 
tal heat capacity of a NS star originates in the nucleons in the 
core. For degenerate fermions of type /, the specific heat per 
unit volume in terms of the dimensionless Fermi momentum 
Xf i — hkf i/niiC is 



1/2 



Then, the contribution of relativistic electrons is 



5.410' 



no 



2/3 



Ti) erg cm ' 



while for non-relativistic nucleons N - n,p is 
c,,j, - 1.6 102» i T.li"- erg cm-^K"' 



mn \ no 



(37) 



(38) 



(39) 



where no - 0.16 fm"^. We include the effect of superfluidity 
through the factor K"' (Levenfish & Yakovlev 1994), which de- 
pends on the pairing state of the nucleons involved ('5o or ^P2)- 
The electron contribution, or that of muons, if present, inside the 
core is, in principle, much smaller, but it dominantes when all 
nucleon species undergo a phase transition to a superfluid state 
(see Fig.|6]l. 

In our model we include the crustal specific heat, which has 
contributions from the neutron gas, the degenerate electron gas 
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3e+18 4e+I8 
V(r)(cm^) 

Fig. 6. Specific heat as a function of the enclosed volume V{r) 
within a radius r, and for different fixed T. We show the con- 
tribution of neutrons (long-dashed line), protons (short-dashed 
line), electrons (dotted-dashed line), and ions (brown solid line). 
Total Ci, is plotted using a black solid line. The dotted lines show 
the non-superfluid case. 



and the nuclear lattice ( Ivan Ripedll991h : it is however negligible 
in comparison to the core contribution, due to the small volume 
of the crust. 

4.4. Neutrino emissivities 

During the first a: 10^ yr, the so-called neutrino cooling era, 
the evolution of a NS is governed by the emission of neutrinos. 
Thereafter, photons radiated from the surface control the evolu- 
tion in the photon cooling era. The path of a NS in a temperature- 
age diagram and the duration of the neutrino cooling era is de- 
termined by the efficiency of the neutrino processes in their in- 
terior. Typically, neutrino emissivities, at high densities, depend 
on temperature to the power of a large number, a nd we akly on 
density; a review is provided bv lYakovlev et al.l (1200 lb . In the 
so-called standard cooling scenario, the total emissivity is dom- 
inated by slow processes in the core, such as modified Urea 
(MUrca) and nucleon-nucleon (N-N) Bremsstrahlung. The min- 
imal cooling model, in which pairing betw een nucleons an d 
the effects of superfluidity are both included dPage et al.ll2004l) . 
is more realistic. If fast neutrino processes, i.e. direct Urea 
(DUrca), occur, the evolution of a NS alters significantly, lead- 
ing to the enhanced cooling scenario. Nevertheless, DUrca only 
operates above a critical proton fraction Yp ^ 0.11, that is 
only reached at high density (4-6) po in the inner core of high 
mass NSs. Since we assume a superfluid core with '5o paring 
for protons and 2 pairing for neutrons, we account for the 
exponential suppression of these processes through reduction 
functions H ( iYakovlev et alj|200 r). We include the Cooper Pair 
Breaking and Formation emissivity (CPBF), although the effec- 
tivity of t his process was questioned both, by observational ar- 
guments (ICumming et al.ll2b06i) and by theoretical calculations 
dLeinson & Perezll2006l) . In the latter work, the authors showed 
that the neutron CPBF emissivity is suppressed, which is 
relevant in the crust, but the proton 'S'o and neutron chan- 
nels, which are relevant in the core, are not seriously altered. 
Including or not this suppression has an effect on the early relax- 



ation of the crust, but has little imprint on the long term cooling 
evolution. 

In our calculations, we consider all relevant neutrino emis- 
sion processes listed in Table [3] which indicates the density and 
temperature dependence of the emissivity for the different pro- 
cesses. The factors that account for further corrections, due to for 
example effective masses and correlation effects, can be found in 
the references listed in Table [3] The table also includes the criti- 
cal proton fraction that is required before some processes can 
operate. In the enhanced cooling scenario, we include the fast 
DUrca process. The efficiency of this fast reaction is exponen- 
tially reduced when superfluidity is taken into account. 



Table 3. Neutrino processes and their emissivities Q in the core 
and in the crust. Third column shows the onset for some pro- 
cesses to operate (critical proton fraction Y^). Detailed functions 
and precise factors can be found in the references (last column). 



Process 



(2[ergcm '] 



Onset 



Ref 



Processes in the core 



MUrca («-branch) 
nn — > pneVe 
pne — > nnVe 
MUrca (p-branch) 
np — > ppeve 
ppe — » npVe 


8 X lO^i K"" n]P Tl 
8 X 10^' n-^" n]l^ Tl 


Yl = 0.01 


1 
1 


NN-Bremsstrahlung 
nn — > nnvv 
np — > npvv 
pp — > ppvv 


7x W^'R""nJ^Tl 
1 X lOF'^ H"" nj^ Tl 
7x WHPi'nJ^Tl 




1 
1 
1 


e-p Bremsstrahlung 
ep — > epvv 


2x\0" n^^^Tl 




2 


DUrca 

n — > peVf, pe — > nv^ 
n -> p^v^,pfi -» nv^ 


AxWH"" n'J^Tl 
AxWn"" n'J^Tl 


y; = o.ii 
y,^ = 0.14 


3 
3 


Processes in the crust 


Pair annihilation 

ee^ — > vv 


9x 1020Fpa„(«<„«,+ ) 




4 


Plasmon decay 
e — * evv 


lxl02»/p,(r,y,) 




5 


e-A Bremsstrahlung 
e(A,Z) -» e(A,Z)vv 


3 X 10'" Lm Zpo n, 




6 


A'-A'-Bremsstrahlung 
nn — > nnvv 


7x lO''^'R""fynlPT^ 




1 


Processes in the core and in the crust 


CPBF 

B + vv 


1 X 10^' n'J' Fa,b Tl 




7 


Neutrino synchrotron 
e — > (B) — > evv 






8 



Ref. (1) 'Yakovle v & LevenfishI |)1995'); (2) 'Maxwell' ('19791); 
(3) E^timer et al. dl991h : (4) Kaminker & Yakovlev ( 199§- 
(5) Yakovlev et al. (2001); (6) Haensel et al. (1996); Kamink er et al] 
a999.) ; (7) Yakovlev et al., a999.) : (8) ,Bezchastnov et al., a997i) 



The neutrino energy losses from processes that occur inside 
the crust are very important at the beginning of thermal evolu- 
tion, during the relaxation stage prior to the core-crust thermal 
coupling. This stage lasts about 10 - 10^ yr and was studied 
in detail by Gnedin et al. (2001). These reactions occurs in a 
wide range of matter compositions, which includes a strongly- 
coupled plasma of nuclei and electrons in the outer layers, a lat- 
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tice of neutron-rich nuclei, up to the crust-core interface with 
abundant free neutrons. Thus, the resulting emissivities are com- 
plicated functions of the temperature and matter composition. 
Considering that the free neutrons in the crust are very likely to 
pair in the '5o state, we account for, in addition, the superfluid 
corrections and the CPFB process in the crust. 

Finally, we regard relativistic electrons that can emit neu- 
trino pairs under the presence of a strong magnetic field, which 
is analogous to the synchrotron emission of photons, because 
our primary goal is to describe the cooling of magnetized NSs. 
This neut rino synchrotron emissiv ity is proportional to the field 
strength (iBezchastnov et all 1 19971) and becomes important for 
B > 10"^ G. 

The emissivities of the most relevant core and crust neutrino 
processes for the minimal cooling scenario are plotted in Fig.|7] 
for three fixed temperatures of 3 x 1 0^ K, 5 x 1 0^ K, and 1 x 1 0^ K. 




le+10 le+U le+12 le+13 le+14 
P (g cm" ) 

Fig. 7. Neutrino processes in the crust and in the core for the min- 
imal cooling for diff'erent fixed T. Lines denote: MUrca (thick 
solid), n-n Bremsstrahlung (long dashed), n-p Bremsstrahlung 
(short dashed dotted), p-p Bremsstrahlung (double dashed dot- 
ted), e-A Bremsstrahlung (long dashed dotted), Plasmon decay 
(solid with cross symbols), CPBF (thin solid) and v-Synchrotron 
(short dashed), assuming a constant field of B = lO'* G. 

At r = 3 X 10^ K, in the early evolution, the plasmon decay 
dominates the neutrino emission in the crust and the MUrca is 
the strongest energy loss mechanism in the core, as seen in the 
upper panl of Fig. |7] At this temperature, neutron superfluidity 
already exists in the crust, and CPBF becomes important near 
the crust-core interface. On the other hand, protons and neutrons 
in the core have not yet condensed into a paired state in a signif- 
icant volume in the core. Later, at intermediate temperatures of 
r = 5 X 10** K (middle panel), plasmon decay is dominant only in 
the outer crust, while electron-nuclei Bremsstrahlung becomes 
more eflicient in a large part of the crust volume. In addition, 
there is an enhancement of the emissivities due to the CPBF, 
at densities between lO'^-lO'* g cm"^. In the core, '5o proton 
superconductivity and neutron superfluidity have a twofold 
effect; suppression of the otherwise dominant processes (MUrca 
and A^-A^ Bremsstrahlung), and enhanced emissivity from CPBF. 
At later stages, when the temperature has fallen to T = 10** K, 



neutrino synchrotron overcomes the other emissivities if a mag- 
netic field of the order of B ^ lO''* G is present. A narrow density 
window is still controlled by CPBF of neutrons in the crust. 



5. Cooling of weakly magnetized neutron stars 

Our discussion in based on two baseline models (see Table 14. il l 
that correspond to the minimal and enhanced cooling scenarios. 
The first case corresponds to low mass NSs in which the central 
density is below the critical density for the onset of the DUrca 
process, which is 2.6 x 10''' g cm"^ for our EoS. The second case 
describes the thermal evolution of a high mass star for which the 
DUrca process operates in a finite volume in the core. For both 
models, we solved the thermal diffusion Eq. (|22] | using several 
magnetic field configurations described in Sect.|2]and the micro- 
physics inputs presented in Sect. |4] including the effects of su- 
perfluidity. We use a two-dimensional numerical grid containing 
350 radial and 60 angular points. 



--^-^ 5.1. Crust formation 



We address the timescale for both the crust formation and the 
growth of the core region where protons are in a superconduct- 
ing state, because the temperature of a NS falls below 10 "'K a 
few minutes after birth. The comparison of these two timescales 
is relevant to understand whether or not there is enough time 
to expel magnetic flux from the core before the crust is formed 
and the magnetic field is frozen into the solid lattice. If this is 
the case, after the crust is formed, the problem can be treated 
by assuming that the magnetic field evolves independently in 
the crust, without penetrating the core, while the thermal evolu- 
tion of the core and the crust become coupled. In contrast, if a 
substantial part of the magnetic flux remains within the core, it 
is probably organized into superconducting flux tubes that have 
a complex interaction with the normal phase. This would be a 
much more difficult problem to solve and the evolution would 
depend on the interaction between the flux tubes and vortices 
and how they become attached to the lattice. The study of such 
a scenario is beyond the scope of this paper, and, for simplicity, 
we assume that either the magnetic field has been completely 
expelled from the core or that it penetrates into the core without 
considering superconducting effects. 

We followed two indicators of the growth of the crust and 
the superconducting core: 

i) the Coulomb parameter, that describes the physical state 
of the ions, defined as 



(Zef 0.23 



kTa; 



(40) 



where o/ = (3/47rn,y is the ion-sphere radius, n/ is the ion 
number density, and T(, is the temperature in units of 10'' K. 
When r < 1, the ions form a Boltzmann gas, when 1 < F < 175 
their state is a coupled Coulomb liquid, and when F > 175 the 
liquid freezes into a Coulomb lattice. The melting temperature 
(Tm) for a body-centered cubic (bcc) lattice corresponds to the 
value at which F = 175. For p - lO''* g cm"-*, we have that 
^ 3 X 10'° K , and the inner crust begins to form at very 
early stages of evolution. We show the evolution of F for the LM 
model in the right panel of Fig. |8] where each line corresponds 
to a different time. The inner crust, up to a radius of 12.4 km, 
has formed completely on a timescale of several hours to a few 
days. To form the outer crust, however, takes much longer, about 
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1-100 yr. The solidification depends, in principle, on the mat- 
ter composition, but we obtained similar qualitative results after 
varying the EoS. 
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Fig. 8. Left panel: Growth of the proton superconducting region 
in the core. T/r^p as a function of r at fixed evolution times. 
Solid (dashed) lines correspond to high (low) T™™ for p '^o, i.e. 
case e (/). Right panel: Crust formation. F vs r at fixed evolution 
times. In both pannels the LM model (minimal cooling) is used. 



ii) the dimensionless parameter T ITc,p for ^Sq proton pair- 
ing; its evolution describes the growth of the superconducting 
region in the core, since for T < T^- j, protons become superfluid. 
We compare two pairing models taken from Table|2] case e with 
high r™" (^ 7 X lO** K) and case / with low T™" (^ 2 x lO** K). 
The evolution of T jT^ p is shown in the left panel of Fig. |8] We 
found that a large part of the core becomes superconducting on 
a timescale that varies from several days to months, depending 
on the pairing model. 

Consequently, we found similar timescales for the formation 
of the solid crust and for the growth of the superconducting core. 
Although it would be interesting to investigate how these two 
processes compete, it is beyond the scope of this work. We as- 
sume that our initial configuration is a NS with a magnetic field 
that remains fixed after the first few days, which is much shorter 
than the overall cooling evolution time. 



5.2. Minimal and Enhanced cooling 

To evaluate whether all microphysics inputs are implemented 
properly in our two-dimensional code, we revisit cooling curves 
for weakly magnetized neutron stars, i.e. with field strengths 
B < 10'^ G. We compare our results considering that, for weak 
fields, the temperature profiles are almost spherically symmetric, 
with previous one-dimensional calculations performed by other 
authors. The most important deviations between models arise, as 
expected, from the underlying microphysics. Major differences 
depend on the occurrence of superfluidity and whether slow or 
fast neutrino emission processes are taking place. We summarize 
these effects below. 

We plot the surface temperature for the minimal and en- 
hanced cooling scenarios in Fig. |9l In both cases, we explored 
different superfluidity models. The major uncertainties come 
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Fig. 9. Cooling curves of weakly magnetized NS with B < 10'^ 
G. Surface temperature (T,) vs age (f) for LM and HM stars. 
n '^0 and p ^Sq are fixed to case a and e, respectively. Solid 
lines show the case without n^P2 gap; for dashed lines n^Pa is 
fixed to case k and for dotted dashed lines to case m. Dotted lines 
represent no superfluidity. 



from the n 2 pairing gap in the core; furthermore, this gap 
is expected to ha ve the strongest impact on the luminosities 
(iPage etal.ll2004 . Hereafter, we fix two models of the n ^Sq 
superfluidity in the crust and the p ^Sq in the core to the cases 
a and e, respectively. We checked that replacing them with the 
other options listed in Table |2] produces slight deviations in the 
cooling curves. We only vary the gap model of the n state be- 
tween three diff'erent limit cases: no pairing (no SF), case h for 
high r™" (^ 6 X lO'' K), and case m for low T™" (^ 2 x 10** K). 

In the early stages of evolution (up to ^ 10^ yr) during the 
initial thermal relaxation of the crust, the main effect of superflu- 
idity is the suppression of the specific heat of free neutrons in the 
crust (see Fig. |6]l, which leads to a faster temperature decrease 
compared to the nonsuperfluid case (dotted lines). The following 
epoch (up to ^ 10^-10^ yr) is controlled by neutrino emission 
from the core. The MUrca and Bremsstrahlung processes (or the 
DUrca process for model B) are exponentially suppressed, in ad- 
dition to the heat capacities of neutrons and protons in the core. 
Nevertheless, core CPFB is important and acts in the opposite 
direction, increasing the emissivities but inside a narrow density 
window. The overall effect is a faster cooling of the LM star. The 
opposite effect is found for the HM star, where a high T™" pair- 
ing of the n ^P2, which covers all the core density region (case 
h), produces a significantly higher Tj with respect to the non- 
superfluid case. If the pairing has a low T™" (case m) or does 
not occupy the whole core volume, then the DUrca process is 
as efficient as in the nonsuperfluid case, leading to a very rapid 
cooling. In the later cooling phase (from 10^-10^ yr), when pho- 
ton luminosity gradually overtakes the neutrino luminosity, the 
most important effect of superfluidity is the reduction of the core 
specific heat that makes the star cool faster. 

In brief, we confirm all previous results for the cooling 
of non-magnetized neutron stars an d do not find any qualita- 
tive difference f rom earlier works (lYakovlev &Pediicki 120041: 
iPageet al.ll2006l) . 
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6. Cooling of strongly magnetized neutron stars 

We study now magnetized NSs with B > 10^^ G. After analyzing 
the effect of superfluidity on the cooling curves, we restrict fur- 
ther study of magnetized neutron stars to two different limiting 
scenarios with fixed superfluidity models, which can be summa- 
rized as follows: 

1. Model A (minimal cooling): a LM star with n^So (case a) in 
the crust, p^So p (case e), and n 2 (case h) in the core. 

2. Model B {enhanced cooling): a HM star with n '^o (case a) 
in the crust, p ^Sq (case e), and non-superfluid neutrons in 
the core 

In this section, the magnetic field given provided by the ini- 
tial model is kept fixed throughout the entire evolution. The ef- 
fect of field decay is separately discussed in the next section. 
As outlined before, one of the most relevant eff'ects of the mag- 
netic field is to reduce the electron thermal conductivity across 
magnetic field lines. Therefore, heat is essentially forced to flow 
along magnetic lines, which results in anisotropic temperature 
distributions. Another important eff'ect is that, as a consequence 
of the different thermal conductivities in the crust and the core, 
their thermal evolution is not always coupled. This effect is also 
magnified by the presence of strong fields. After the initial fast 
transient in which large gradients are allowed due to the reduced 
thermal conductivity at high temperature, there are different pos- 
sible evolutions depending on the magnetic field geometry. Since 
the field lines close to the poles are essentially radial, in general 
the magnetic poles are thermally connected with the core and 
reflect its temperature. In contrast, the equator is insulated by a 
magnetically-induced thermal wall due to large tangential com- 
ponents. Its evolution is thus almost independent of that of the 
core. We discuss first our results for purely poloidal fields and 
then consider the effect of toroidal fields. 



6.1. Purely poloidal magnetic fields 

In Fig. [TOl we plot temperature profiles across the star, for the 
Model A, as a function of the density, and for different evolution 
times. The magnetic field is confined to the crust (PC; and Bg 
as in Fig. [U and B = 5 x 10'^ G. In the early stages, the pole 
cools down in a similar way to the core and its temperature is 
almost identical to that of the non-magnetized case. On the other 
hand, the equatorial region is decoupled and shows a different 
evolution. Since the crustal heat cannot be released inwards, into 
the core, where neutrino emission is an efficient cooling mecha- 
nism, it remains warmer for a longer time, typically few 10^ yr 
(Fig. [TOl left panel). At intermediate ages, a nearly isothermal 
state is reached and the crust and the core evolve together, ap- 
proximately from 10^ yr to lO'* yr (Fig.fTO] right panel). At the 
late evolution (10^-10^ yr), photon emission from the surface is 
the most efficient way to radiate energy and the initial situation 
is reversed: since the equator cannot be refed by the relatively 
warmer core, it becomes cooler than the pole. We obtained the 
same qualitative results for Model B. 

In Fig. [TTl we show the magnetically induced anisotropic 
temperature distribution for the same model. The upper panels 
are the usual cooling curves (temperature vs. time), which dis- 
play the evolution of T/, at the magnetic pole and the equator, for 
two different field configurations: core dipolar (CD) and poloidal 
confined to the crust (PC). For the latter configuration we show 
results for two field strengths. The cooling curves do not show 
large deviations between models, although the difference with 
respect to the non-magnetized case becomes larger with increas- 
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Fig. 10. Evolution of the temperature as a function of the density 
of a LM magnetized NS with B = 5 x lO'"' G (PC configuration). 
T at the pole is shown with thick solid lines and at the equator 
with thin solid lines. Evolution times are indicated near the lines. 
The B - Q case is plotted with dashed Unes. 



ing magnetic field strength. The lower panel shows the corre- 
sponding angular distribution of Ti,, normalized to its value at 
the pole, for three different ages. At f » 500 years, we find 
that, for all models, the crustal temperature at the pole is smaller 
than at the equatorial region. We refer to an inverted tempera- 
ture distribution, in such cases, when we find cooler polar caps 
with a warmer equatorial belt. The occurrence of this inverted 
profile is model independent while its duration and the degree of 
anisotropy reached depend on the details of the magnetic field 
geometry and strength. The equivalent results for Model B are 
shown in Fig. [12] For example, for model B, in which the DU 
process is allowed, the star cools faster and its interior reaches 
higher values of the magnetization parameter, making the in- 
verted temperature profile more pronounced {T^^ ^ 27'^°''^). For 
crustal magnetic fields the inverted profile can be maintained for 
longer times (^ 10^ yr) than for the core dipolar case that be- 
comes isotropic at about 10^ yr (lower panel). During late evo- 
lution, at 10^ yrs, the usual temperature distribution is found: a 
hot polar cap with a cooler equatorial belt. 

We reiterate that Tt, is the temperature at the bottom of 
the envelope, corresponding to our outer point of integration 
at p a: 10^ g cm"^, and the blanketing effect of the envelope 
should be taken into account before comparison with observa- 
tions. To translate the temperature at the base of the envelope 
to the surface temperature Tj, we assumed a magnetized enve- 
lope as described in Sect. 13.21 taking into account the angle that 
the magnetic field forms with respect to the normal to the sur- 
face. In Fig. [T3] we plot Tj and T JT^"^'' vs. age for the same 
three cases as in Fig. [TT] We note that the anisotropy found at 
the level of Ti, does not automatically produce a similar surface 
temperature distribution: the blanketing effect of the envelope 
overrides the inverted temperature distribution found at interme- 
diate ages. We see that the equator remains always cooler than 
the pole, and only at early times and for strong fields do we find 
larger surface temperatures in middle latitude regions. Another 
important result is that the degree of temperature anisotropy at 
the level of the crust is always rather small for magnetic fields 
penetrating into the core (i.e. CD), which causes a similar sur- 
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Fig. 11. Cooling of strongly magnetized NSs for the Model A. 
Upper panel: Ti, vs t, at the pole (thick solid lines) and at the 
equator (thin solid lines). Two field configurations are shown: 
CD (left, for B = 5 X 10'^ G) and PC (middle, for B = 5 x 10'^ G, 
and right for B = 5 x lO''* G). Dashed Hnes indicate the B = 
case. Lower panel: TblT'^"^^ vs. the azimutal angle 9 for three 
fixed evolution times 20 yr (solid lines), 500 yr (dashed lines), 
and 10^ yr (dotted dashed lines). Similar field configurations as 
in the upper panel are shown. 



face temperature distribution at all times and independently of 
the magnetic field strength. Crustal confined fields, in contrast, 
allow for non-uniform temperatures at the base of the envelope, 
leading to temporal variations in the surface temperature distri- 
bution during the NS life. Since we are interested in the models 
with the largest variation of temperature, hereafter we only con- 
sider crustal confined fields. 

If we c ompare our temperature angular distribution to for- 
mer results jGeppert et al.ll2004t iPerez-Azorm et al]|2006ah . we 
find that our late time profiles coincide qualitatively with the 
stationary solutions obtained in previous works. However, the 
temperature distributions at early times are quite different. The 
reason is that stationary solutions cannot describe properly the 
temperature distribution in young NSs, because a NS is evolving 
and changing its thermodynamical conditions faster than, or on 
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Fig. 12. Same as Fig.[n]but for the Model B. 



a similar timescale, to the time needed to reach the stationary 
state. 

6.2. Effect of toroidal fields 

Despite our lack of direct information about the magnetic field 
geometry inside a neutron star, there is some agreement that sev- 
eral independent mechanisms can create strong toroidal fields, 
such as differential rotation during core collapse (' Wheeler et al.l 
^002), or proto-neutron star dynamo. Hence, it is natural to 
investigate the effect of toroidal components on the surface 
temperature distribution. In particular, the inclusion of toroidal 
components was used to explain t he small hot emitting ar- 
eas observed in som e isolated NSs ( Perez- Azorm et ani2006bt 
iGeppert et al.l |200"6|) . These works concluded that the surface 
temperature is determined more by the geometry rather than by 
the magnetic field strength. With this motivation, we include a 
toroidal component in our models and study its influence on the 
results. 

In the remainder of this section, we focus on the effect of the 
toroidal component on model A, but our qualitative conclusions 
can be generally extended to model B. In Fig. [141 we compare 
the cooling curves in the upper panel, and the angular tempera- 
ture distribution in the lower panel, for B = 5 x lO''* G, for the 
different toroidal field configurations. The first conclusion is that 



14 



Deborah N. Aguilera et al.: 2D Cooling of Magnetized Neutron Stars 




log t (yr) log t (yr) log t (yr) 



'0.6- 



0.2 



: B = Sxlo'^G (CD) 


: B = 5xio''G(PC) : 


: B = sxio'^G (PC) : 




500 yr I 


- 500 yr 




V \ / -'^ " 

\ \ / / _ 
\ / '' 
^ \ \ f 

: 10 yA\ // : 






Ml' 

M /' - 
H ' 


s * \ / '' 
10 yA / - 




W I 


^ 1 '' - 
r 'j r 


' < < . 1 < < < 1 < . . 


" , , , 1 , , , 1 , , , r 


" , , , 1 , , , 1 , , , r 



1 2 
e (rad) 



1 2 
e (rad) 



1 2 3 
e (rad) 



Fig. 13. Cooling of strongly magnetized NSs for the Model A. 
Same as Fig.[TT]but for the surface temperature T s (upper panel) 
and for Ts/T'^°^'^ (lower panel). 



the presence of crustal confined toroidal fields (TCI and TC2) 
does not significantly change the results obtained with purely 
poloidal fields (dotted lines in Fig.[T4ll. We omitted model TC2 
in the upper panel because it was indistinguishable from TCI; 
minor differences are visible only during the first 100 years of 
evolution (see lower panel). We found that the FF configuration 
exhibits larger Tb than the other models in the late evolution 
(f > 10^ yr) because the heat transport was suppressed by the 
toroidal component extended through the envelope, and the in- 
sulating effect was more pronounced. The reason for the large 
differences between models with toroidal magnetic fields con- 
fined to the crust and the FF model, was the difference in the an- 
gle that the magnetic field forms with the normal to the surface. 
According to Eq. dSTT i. the more tangential the field (the larger 
(fi), the smaller Ts for a given Th- In general, we expect that the 
presence of toroidal fields extended to the envelope and magne- 
tosphere results in lower surface temperatures in the equatorial 
region. During the neutrino cooling era, the polar region remains 
as hot as in the poloidal case, but during the photon era, due to 
the reduced photon luminosity, the star cools more slowly and 
the pole remains warmer 

As we can see in the lower panel of Fig. [14] for the FF case, 
the toroidal field maintains, during the entire evolution, a cooler 
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Fig. 14. Cooling of strongly magnetized NSs with toroidal fields. 
Upper panel: Ti, vs t (on the left) and Ts vs t (on the right) at 
the pole (thick lines) and at the equator (thin lines). Three field 
configurations with B = 5 x 10'^ G are shown: PC (dotted lines), 
FF (solid lines) and TCI (dashed lines). Lower panel: Ts/Ts°^'' 
vs 6 for three fixed evolution times. Same field configurations as 
in the upper panel are shown. 



and more extended equatorial belt, while the hot polar region is 
shrunk in comparison to the other models. Defining the angular 
size of our polar cap by the angle at which the radial compo- 
nent of B becomes larger than the tangential component, this is 
> Bg, for a purely poloidal configuration, which implies a 
hot area of about ^ 40° - 60°. For a FF model, this condition is 
reached when B^. > Bg + B^, which provides a smaller angular 
size of about 10°, which agrees with the estimated emittin g area 
for some isolated neutron stars (^ Perez-Azorm et al ]|2006b). The 
same comment made at the end of the previous subsection about 
the comparison with stationary models is valid when comparing 
these r esults to stationary temperature distributions wi th toroidal 
fields jPerez-Azorm et al.ll2006aHGeppert et al.ll2006l) . 

7. Magnetic field decay 

In the previous section, we discussed the impact of strong mag- 
netic fields on the cooling and the temperature distribution of 
NSs, keeping the field strength and geometry fixed throughout 
the entire evolution. But the existence of crustal confined fields 
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supported by crustal currents is inconsistent with the assump- 
tion of non-evolving fields. Currents in the crust dissipate in a 
relatively short timescale, which may vary depending on the in- 
teraction of electrons with the lattice in different crustal regions. 
In any case, this leads to dissipation of the magnetic field on 
timescales at least comparable, if not shorter, than the cooling 
timescale. This effect is important while the crust is still hot be- 
cause of the large temperature dependence of the electrical re- 
sistivity. At late times, that is after a few 10* years, when the 
crust temperature drops below 10' K, the conductivity increases 
significantly, although limited by electron-impurity or phonon- 
impurity scattering, and the magnetic field decays on a much 
longer timescale. 

Since we are interested mostly in the evolution of NSs while 
their temperature is sufficiently high for them to be visible as 
thermal emitters, the effect of Joule heating by magnetic field 
decay cannot be ignored. 
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7. 1 . Effect of Joule heating on the cooling of magnetized NSs 

Based on more det ailed works stu dying the magnetic field evolu- 
tion in NS's crusts (iPons & Gepp ert 2007), we assumed the sim- 
plified form of field decay provided by Eq. (fTTT i. and we chose 
the model TCI as representative of the type of fields expected to 
arise from those simulations. Although we center our discussion 
on the particular case of model A with TCI, we stress that our 
results depend qualitatively neither on the particular NS model 
(equation of state, crust size, etc.) nor the choice of radial depen- 
dence of the toroidal component. 

Having fixed our background NS model and field geometry, 
we varied the parameters that describe the typical timescales for 
Ohmic dissipation and a fast initial decay induced by the Hall 
drift. In Fig. [15] we show the cooling curves for three different 
pairs of values (rohm, THaii) = {(10^ 10^); (10^ 10^); (10^ 10^)} 
yr, represented by solid lines, dashed lines, and dash-dotted 
lines, respectively. For comparison, the dotted lines show the 
evolution with constant field for the same initial field {Bq - 
SxlO'^G). 

The decay of such a large field has an enourmous effect upon 
the surface temperature; due to the heat released, the temperature 
remains far higher than for a non-decaying magnetic field. The 
strong imprint of the field decay is evident for all pairs of param- 
eters chosen. We note that the temperature of the initial plateau 
is higher for shorter THaii, but the duration of this stage, which 
has almost constant temperature, is also shorter This is a conse- 
quence of releasing a similar amount of heat in a shorter time: at 
t = THaii, B has decayed to about half of its initial value and three 
quarters of the initial magnetic energy has dissipated. By reduc- 
ing THaii we can therefore maintain higher temperatures, but for 
shorter times. After f = THaii, there is a noticeable drop in Ts due 
to the transition from the fast Hall decay to the slower Ohmic 
decay. 

The insulating effect of tangential magnetic fields operates 
in both directions: in the absence of additional heating sources, 
it decouples low latitude regions from the hotter core resulting 
in lower temperatures at the base of the envelope; conversely, if 
heat is released in the crust, it prevents extra heat flowing into 
the inner crust or the core where it is more easily lost in the form 
of neutrinos. Indeed, our simulations that include Joule heating 
systematically indicate the presence of a hot equatorial belt at 
the crust-envelope interface. Kaminker et al. (2006) studied the 
effect of a localized heat source at different depths inside a NS. 
They concluded that only a heat source very close to the stellar 
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Fig. 15. Cooling of strongly magnetized NSs with Joule heat- 
ing with Bo = 5 X 10'"^ G. Upper panel: vs f at the pole 
(thick lines) and at the equator (thin lines) for three pairs 
of values (rohm, THaii): (10^ lO"*) yr (soHd Hnes), (10^ 10"*) yr 
(dashed lines), and (10^, 10^) yr (dotted dashed lines), respec- 
tively. Lower panel: Ts/T'^°^^ vs for three fixed evolution times. 



surface can have observational consequences. In this work, we 
find evidence for a far more important effect on the surface tem- 
perature. The main reason for this apparent discrepancy is that 
our cooling models are two-dimensional and include the insu- 
lating effect of the strong tangential field in the c rust, as opposed 
to the one-dimensional simulations studied by iKaminker et al.l 
(I2006h . However, as discussed in the previous section, this in- 
verted temperature distribution at the level of the crust is not 
necessarily visible in the surface temperature distribution be- 
cause it is filtered by the magnetized envelope. An analysis of 
the angular temperature distribution shown in the lower panel of 
Fig. [Ts] shows an interesting feature related to the heat deposi- 
tion: the development of a middle latitude region hotter than the 
pole at relatively late stages in the evolution (f ^ 10'*, 10^ yr). 
For a wide range of parameters we found this hot area. It would 
have implications for the light curves of rotating NSs, that will 
differ substantially from the light curves obtained with a typi- 
cal model consisting of a hot polar cap with a cooler equatorial 
region. 
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7.2. The hidden direct Urea process ? 

We conclude our discourse on the important impact of magnetic 
field decay in the cooling history of a neutron star by reconsid- 
ering the enhanced cooling scenario, in which the DUrca effi- 
ciently cools the star very quickly. In Fig. [16] we compare our 
results for the cooling of low and high mass NS (Models A and 
B), with and without magnetic fields. Neglecting the effect of 
magnetic fields, the differences between the fast and slow cool- 
ing scenarios (short dashed and dotted lines, respectively) are 
clearly evident, although they can be reduced by strong super- 
fluidity. We consider, however, a limiting case that experiments 
a rapid cooling, which has no superfluidity in the inner core, to 
observe the significance of the magnetic field effects on the sur- 
face temperature (thick sohd and dashed lines). As we can see 



pole 




log t (yr) 

Fig. 16. Comparison of the fast and slow cooling including Joule 
heating. vs age for Model A (solid lines) and for Model B 
(dashed lines). The temperature at the pole (thick lines) and at 
the equator (thin lines) is shown. The initial field is Bq = 5 x 
IQ^'^G and the field decay rate is of rohm = 10^, THaii - 10"* yr. 
The B = case is shown in thin long (short) dashed lines for 
Model A (B). 

in Fig. [16] in a NS born with a field of 5 x 10'"* G that decays 
about one order of magnitude during the first million years of its 
life, it is hard to distinguish whether or not a fast neutrino emis- 
sion process is active. The surface temperature, both at the pole 
and the equator, is essentially determined by the magnetic field 
geometry, strength and decay rate. Only at late times, the dif- 
ferences between models with and without DUrca become vis- 
ible, but still the variations between models with different field 
strengths can be larger than the differences stemming from the 
fast neutrino cooling process. This interesting result could im- 
ply that we need to reconsider the observations because a fast 
neutrino cooling process may well be triggered inside neutron 
stars but hidden by the magnetic field. A detailed analysis of 
the different possibilities of fast cooling (hyperons, quarks, pure 
nucleonic matter with large symmetry energy) and comparison 
with the observations is beyond the scope of this paper, but the 
problem is thought-provoking. Our first results indicate that di- 
rect URCA may be veiled by magnetic fields and this scenario 
may not be properly identified as fast cooling (Aguilera, Pons & 
Miralles, 2007, in preparation). 



8. Conclusions 

We have presented a thorough study of the thermal evolution 
of neutron stars including some of the most intriguing effects 
of magnetic fields. Our results were based on two-dimensional 
cooling simulations of realistic models that account for the 
anisotropy in the thermal conductivity tensor. In the first part of 
the paper, we revisited the classical scenario with low magnetic 
fields and presented the input microphysics, working assump- 
tions, and the baseline models. As an interesting byproduct, we 
reconsidered the growth of the crust and of the superconduct- 
ing region in the NS core, and found that there are situations in 
which both growth rates are comparable. The main body of the 
work was aimed at the discussion of the two principal effects of 
magnetic fields: the anisotropic surface temperature distribution 
and the additional heating by magnetic field decay. We found 
that, even for purely dipolar fields, an inverted temperature dis- 
tribution is plausible at intermediate ages. Thus the surface tem- 
perature distribution of neutron stars with high magnetic fields, 
even in the axisymmetric case, may be quite different from the 
model with two hot polar caps and a cooler equatorial region. 
The irregular light curves of some isolated neutron stars, for in- 
stance RBS1223, (Schwope et al. 2005; Kaplan & van Kerkwiik| 
|2005|) are an indication of such complex structures. 

The main result of this work is that, in NSs born as magne- 
tars. Joule heating has an enormous effect on the thermal evo- 
lution. Moreover, this effect is important for intermediate field 
stars. If the magnetic field is supported by crustal currents, this 
effect can reach a maximum because two combined factors en- 
hance the efficiency of the heating process: i) more heat is re- 
leased into the crust, in the regions of higher resistivity close to 
the surface, and ii) large non radial components of the field chan- 
nel the heat towards the surface, instead of being lost by neu- 
trinos in the core. As expected, it becomes clear that magnetic 
fields and Joule heating are playing a key role keeping magne- 
tars warm for a long time, but it is likely that the same effect, al- 
though quantitatively smaller, must be considered in radio-quiet 
isolated NSs or high magnetic field radio-pulsars. 

Another aspect that should be considered when we try to 
explain observations using theoretical cooling curves is that for 
many objects the age is estimated assuming that the loss of an- 
gular momentum is entirely due to dipolar radiation from a mag- 
netic dipole (spin-down age). In the case of a decaying magnetic 
field, the spin down ge e, seriously overestimates the true age 
(iGunn & Ostrikerlll970 ). Therefore, the cooling evolution time 
should be corrected, according to the prescription for magnetic 
field decay, to compare our model accurately with observations. 
A detailed comparison of the cooling curves obtained in this 
work with observational sources can be found in Aguiler a^ et al.l 
(120081) . 

Our last striking remark is that the occurrence of direct 
URCA or, in general, fast neutrino cooling in NS may be hidden 
by a combination of effects due to strong magnetic fields. Our 
conclusion is that the most appropiate candidates to monitor as 
rapid coolers are NSs with fields lower than lO'-' G. Otherwise, 
we may be misled in our interpretation of the temperature -age 
diagrams. 

The main drawback of our work is that we are not yet able 
to return a fully consistent simulation of the magneto-thermal 
coupled evolution of temperature and magnetic field. In the near 
future, we plan to extend this study by coupling our thermal dif- 
fusion code to the consistent evolution of the magnetic field in 
the crust given by the Hall induction equation. That approach 
will permit the accurate evaluation of the heating rates, includ- 
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ing the non-linear effects associated with the Hall-drift in the 
NS crust. We believe, however, that the phenomenological pa- 
rameterization employed in this paper, reproduces qualitatively 
the results expected in a real case. We have provided another step 
towards understanding the cooling of neutron stars, by pointing 
out a number of important features that must be more carefully 
considered in future work. 
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